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Abstract In this paper we investigate the use of Richardson extrapolation to esti- 
mate the convergence rates for numerical solutions to wave propagation problems 
involving discontinuities. For many cases, we find that the computed results do 
not agree with the a-priori estimate of the convergence rate. Furthermore, the esti- 
mated convergence rate is found to depend on the specific details of how Richard- 
son extrapolation was applied, in particular the order of comparisons between 
three approximate solutions can have a significant impact. Modified equations are 
used to analyze the situation. We elucidated, for the first time, the cause of appar- 
ently unpredictable estimated convergence rates from Richardson extrapolation in 
the presence of discontinuities. Furthermore, we ascertain one correct structure of 
Richardson extrapolation that can be used to obtain predictable estimates. We 
demonstrate these results using a number of numerical examples. 
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1 Introduction 

Estimating the error in numerical approximations to solutions of partial differ- 
ential equations is important for many reasons. In order to be useful, numerical 
simulations should be accurate in some measurable norm. A particular application 
could require the absolute error to be less than a certain level. Other applications 
might use an estimate of the numerical error to help guide decisions about the 
most cost effective way to spend scarce resources, for example to choose between 
higher resolution simulations, to include more physical processes into the model, 
or to produce more samples for a statistical analysis. Still other applications may 
use estimates of the error directly for uncertainty quantification purposes. There 
are many approaches to error estimation found in the literature. Intrusive tech- 
niques such as adjoint error estimators [Ul2], error transport [HH], or finite element 
residual and recovery methods [5] are extremely powerful. However, because they 
are intrusive they require access and modification of the source code. This is often 
not possible for theoretical, practical, or sometimes legal reasons. 

Non-intrusive techniques are those that require only the ability to produce 
multiple simulation results, but do not require modifications of the source code. 
Error estimation through Richardson extrapolation is one commonly used non- 
intrusive technique and essentially relies on asymptotic properties of numerical 
approximation. Asymptotically correct in this context refers to the fact that ex- 
pending a certain amount more computational effort yields a predictable increase 
in the accuracy of the result. Richardson extrapolation estimates can be based 
on varying the order of approximation, varying the resolution of the grid, or a 
combination of both. For many applications, Richardson extrapolation has been 
shown to yield very good results. However, the behavior of Richardson extrapola- 
tion error estimates for simulations of solutions with discontinuities is known to be 
problematic [6]. Discontinuous solutions are common in many physical systems, 
e.g. near material interfaces in electromagnetics or shock and contact waves for 
fluid mechanics, and so understanding the behavior or Richardson extrapolation 
in theses cases can be important. There have been attempts to introduce richer 
ansatz when using Richardson extrapolation in the presence of discontinuities, and 
these have yielded varying degrees of success. However, there has been very little 
progress on understanding the fundamental sensitivity of Richardson extrapola- 
tion error estimates in the presence of jumps. The current work is a step toward 
filling that gap. 

In this paper, we investigate one particular realization of Richardson extrapola- 
tion error estimation for linear wave propagation with discontinuous solutions. Lin- 
ear jumps are important in their own rite, for example wave propagation through 
solids [3[8], and for their significant role in nonlinear problems, for example contact 
and slip surfaces in the Euler equations [9, 10 . In this manuscript, we build on the 
previous work of [11] which analyzed convergence rates for approximate solutions 
of linear advection where the exact solution contained jumps. That work used 
modified equations to argue that the expected rate of convergence for a nominally 
p** 1 order method in the presence of a linear jump discontinuity is p/(p+ 1). In the 
current work we use the structure of the modified equation solutions to discuss 
the expected behavior of Richardson extrapolation error estimates. We show that 
under certain conditions one can expect to obtain the p/(p+ 1) rate. In addition, 
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we show why the method can fail to obtain the correct result if these conditions 
are not met. 

The remainder of this paper is structured as follows. Section [2] discusses some 
preliminaries, and provides a very brief overview of the Richardson extrapolation 
technique for error estimation. A motivating example problem demonstrating the 
difficulty associated with using Richardson extrapolation error estimation for lin- 
ear jumps is presented in section [3] Section [4] gives an analysis of a simplified 
model problem of ID linear advection. In Section |4.1| we apply Richardson extrap- 
olation to approximations generated by a first-order upwind method, and discuss 
the results. The technique is found to be effective for this case, and an analysis 
explaining this surprising result is presented. That analysis is extended in Sec- 
tion [472] to discuss the case of high-order linear schemes. This analysis reveals that 
one particular instantiation of Richardson extrapolation produces the expected 
convergence rate while others may not. Section [4.4| demonstrates the theory for 
upwind discretizations of order 2, 4, and 6, as well as the case of a high-resolution 
nonlinear TVD discretization. Additional details of the inner workings are pre- 
sented for the second-order case. In Section[5] the motivating example of Section[3] 
is revisited, and the conclusions from the ID analysis are shown to hold even for 
this more complex 2D case. Concluding remarks are presented in Section [6] 



2 Preliminaries and Richardson extrapolation for smooth problems 

Richardson extrapolation is a commonly used technique for error estimation, and 
many variations exist. For a good overview of the technique refer to [6]. Here we fo- 
cus on one particular approach to Richardson extrapolation which uses numerical 
approximations at three grid resolutions obtained using the same numerical tech- 
nique. Even within this particular flavor of the approach, there are essentially three 
possible realizations. In this section we review the approach and present the three 
choices. In what follows, we consider numerical approximations to the solution of a 
partial differential equation (PDE) . Boundary conditions are an important aspect 
of many numerical simulations, but are not critical to the present discussion. Thus 
consider the spatial domain x £ O, and introduce a spatial discretization with 
uniform grid spacing h. 

Consider a set of numerical approximations given by Uh M (x, t) « M e (x, t) where 
u e (x, t) is the exact solution, and indicates the size of the mesh. We consider 
performing an estimate at some time t = tf, and whenever the time argument is 
not included, it is assumed to imply t = tt (i.e. u(x) = u(x, tf)). Let an estimated 
convergence rate be denoted by H(uh % ,Uh 2 ,Uh 3 ) where the various uf lM are nu- 
merical approximations obtained using grid spacing h^.j, and 7t(uf ll ,Uf l2 ,U} l3 ) = a 
is the solution of the scalar equation f(a;ui ll ,ui l2 ,U) l3 ) = 0, where 

f(a;u hl ,u h2 ,u h3 ) = T — —r- - — — (1) 

For the purposed of the remainder of this paper we will assume that ||-|| indi- 
cates a discrete approximation to the L\ norm. This is the norm which is most 
often considered when discussing hyperbolic equations with discontinuities. For a 
given set of three numerical approximations, there are essentially three distinct 
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ways that the estimate can be computed 1Z(uf ll ,iif l3 ,Uh 2 ), TL (uh x , u h 2 > u h 3 ) > and 
lZ{u h2 ,u hl ,u h:i ). As shown below, this distinction is irrelevant for smooth prob- 
lems. However, it will become important for solutions with discontinuities. For 
smooth problems, the assumption underlying basic Richardson extrapolation er- 
ror estimation is that a given numerical approximation u) lM (x, t) is related to the 
exact solution « e (x, t) as 

u h M {*-,t) = UeCM) + c(x,t)h P M + . . . 

where p is the formal order of accuracy of the approximation, and c(x, t) is an order 
one function which is independent of the mesh parameters. High-order terms are 
ignored, and the difference between approximations at two resolutions is 

u/n(x) - Uh 2 (x) = M e (x) + c(x)h p - lte(x) - c(x.)h P 
= c(x) (h\ - h p ) . 

For any three resolutions then we find that 

K t (x) -~"ft 2 (x) 

IK 2 ( x ) - u ^ 3 ( x ) 



Assuming hi ^= h% ^ /13, Equation ([2| and its counterparts can easily be used to 
show that 71 , u, l3 ,u h2 ) = Tl(u hl , u }l2 , u, l3 ) = TZ {u h . 2 ,u hl ,u ha ) = p. Note that in 
principle one can then use the computed convergence rate a in order to estimate 
the exact solution and obtain field estimates of the error. Such an approach is 
presented in detail in [12] and [13] . Also note that there is the possibility multiple 
roots in ([T]), but this situation is easily recognized in practice and so we do not 
discuss this further. 



3 A motivating example 

In order to motivate the need for detailed understanding of Richardson extrap- 
olation error estimators for linear wave propagation problems involving disconti- 
nuities, we revisit the recently published article [Hj. In that paper, the authors 
develop a new class of numerical methods for second-order wave equations. In two 
space dimensions, schemes of order 1, 2, 4, and a nominally 2 order nonlin- 
early limited scheme are presented. For the linear schemes the order of accuracy is 
proved. In order to demonstrate the properties of the methods, a number of test 
problems are presented. Where the exact solution is known and sufficiently smooth, 
the theoretical accuracy is confirmed. This provides confidence that the computer 
code correctly implements the numerical methods. However, for arguably the most 
interesting test problem, the exact solution is not known and convergence was 
judged only visually. In fact, the reason Richardson extrapolation was not used to 
quantify convergence for that case was the inherent inconsistency produced by the 



i] C (x) (h\-m 

\\c(x)(h%-hl)\ 
l|c(x)|||^-^| 
||c(x)|| \h p 2 -h p 3 \ 

\t£-hS\- 



(2) 
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method, and the subsequent difficulty in interpreting the results. That difficulty 
is reproduced here. 

Consider the initial boundary value problem 



dt 2 4 dx 2 + dy 2 



(x,y) 6 (-7T,7r) x (-2tt,0), 



du , 



u(x,y,0) = Uo(x,y), ~(x,y,0) = v (x,y), 
du du 

— (x,0,t) = 0, — (x, ~2n,t) = 0, u(x - 7r, y, t) = u(x + it, y, t) 
dy dy 



(3) 
(4) 
(5) 



where 



vo(x,y) = 0, u (x,y) 



1 if x 2 + {y + 2f < 1, 
otherwise. 




Fig. 1 Time evolution of the two-dimensional top-hat using a fourth-order accurate upwind 
method and 6400 points in each coordinate direction. See (14) for details concerning the method. 



Figure [T] shows the evolution of the numerical solution computed using the 
fourth-order accurate scheme from [T3] at times t = 0, t = 1.5, and t = 3.0. The 
solution was computed using 6400 points in the x- and y-directions. In order to 
judge convergence one could apply Richardson extrapolation as discussed in Sec- 
tion [2] Table [T] presents such a study using uniform refinement and shows results 
for the three variants of Richardson extrapolation as applied to approximations 
from numerical methods of order 1, 2, 4, and a high-resolution nonlinear variant of 
nominally second-order. As discussed in [ll|H4j . the expected convergence rate for 
this case is p/(p+ 1), where p is the nominal order of the method. The challenge 
in interpreting the results in Table [l] is apparent as there is wide variations in the 
approximated convergence rate depending on which ordering within the Richard- 
son extrapolation scheme is used. Indeed the fourth-order scheme yields results 
between —.5 and 3.41, despite the fact that reasonably highly resolved simulations 
are being considered. Furthermore, we note that the results presented in the first 
column are in reasonable agreement with the expectation. 
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scheme 


ft (u hl ,u hl / 2 ,u hl/A ) 


ft(u hl ,u hl /4,u hl / 2 ) 


ft (uh 1 /2,u hl ,u hl/4 ) 


first order 


0.44 


0.40 


0.50 


second order 


0.53 


-0.27 


1.92 


high resolution 


0.62 


0.52 


0.79 


fourth order 


0.64 


-0.50 


3.41 



Table 1 Estimated convergence rates for the top-hat problem using upwind schemes and 
the three variants of Richardson extrapolation. The base resolution uses 400 points in each 
coordinate direction, and uniform refinement is carried out using a ratio of 2. Results are 
presented for numerical methods of order 1, 2, 4, and a high-resolution nonlinear variant of 
nominally second-order. 



4 A model problem with discontinuity 

We now seek to understand the nature of Richardson extrapolation error estima- 
tion for problems with discontinuities, or other self similar behavior, using a simple 
model problem. Consider the one dimensional linear advection equation 

d d 

—u(x,t)+a—u(x,t) = (6) 

with constant advection velocity a > 0. A canonical model problem with disconti- 
nuity can be defined using the initial conditions 

u{x,0) = ■{ (7) 
v ' I ur for x > 0. 

The method of characteristics is used to define the exact solution for all t > as 
u(x,t) = u(x — at,Q), which applies also to discontinuous solution profiles using the 
notion of weak solutions |151l9] . We note here that the expected order or accuracy 
for all numerical methods considered in this paper are confirmed using known 
smooth solutions. 



4.1 First Order Upwind Discretization 

Consider the first-order accurate explicit upwind scheme 

n+l n \ r n n i /o\ 
u i = u i ~ X l u i - u i-l\ (°) 

where w™ is a numerical approximation to u(xi,t n ) and the so called CFL number is 
A = ^t—. The computational domain [a;^, xg\ is a truncation of the infinite domain, 
and has been discretized with Xi = xl + ih where h = (xr — xl)/(N — 1) and N 
an integer. Zero gradient conditions are applied at domain boundaries, but the 
nature of the specific problem being studied makes the details of these artificial 
boundary conditions unimportant. Time has been discretized as t n = nAt with 
initial conditions u(x, 0) being given at t = with = — 1 and ur = 1. Numerical 
stability is obtained for A < 1. 
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4-1.1 Richardson extrapolation error estimation 

One can perform an estimate of the convergence rate using Richardson extrapola- 
tion and simply ignore the fact that the assumptions underlying the approach are 
not strictly valid for cases with discontinuities. We set a = 1, choose a computa- 
tional domain with [a^x/j] = [— 7r,7r], integrate to a time tf = 2, and use A = 0.6. 
Note that in order to ensure that the initial condition is applied consistently for all 
cases, the grid point located at x = is explicitly taken to have the initial value ut,. 
A series of approximations is generated using a uniform refinement process with a 
ratio < r < 1 (i.e. hi = rh% and /13 = r 2 hi), and starts with 51201 points in the 
domain (i.e. h\ = 51 2 2 ^ ). Table [2] shows the results using the three basic variants 



r 


{ u h 1 ,u rhl ,u r 2 hl ^ 




(u r h 1 ,U hl ,U r 2 hl ' S j 


1 

2 


0.50 


0.50 


0.50 


2 
5 


0.50 


0.50 


0.50 


1 

3 


0.50 


0.50 


0.50 


2 
7 


0.50 


0.50 


0.50 


1 

4 


0.50 


0.50 


0.50 



Table 2 Estimated convergence rates for the first-order upwind scheme for a solution with 
a discontinuity. The base resolution uses 51201 points, and uniform refinement is carried out 
using a ratio of r. Results using the three independent variants of Richardson extrapolation 
are presented in the various columns. 



of Richardson extrapolation, and various choices of the uniform refinement ratio 
r. The table makes clear that the estimated convergence rate is 0.5 for any choice, 
which agrees exactly with the expected convergence rate for a first-order scheme 
with a linear jump |llj . 



4-1.2 Explanation of the result 



The fact that Richardson extrapolation seems to work, in terms of convergence 
rate estimates, for the first order upwind scheme even with a discontinuous solution 
is surprising. In order to understand this result we extend the analysis in [11] . The 
approach makes use of modified equation for a more complete understanding of the 
behavior. The modified equation is a continuous PDE whose solution describes the 
approximate behavior of the well resolved components of the discrete solutions, and 
is derived by substituting continuous functions U(x,t) into the discrete equation 
Q by setting uf = U(xi,t n ), and expanding all terms in Taylor series about the 
point (x,t) = (xi,t n ). For the first-order upwind scheme the result is 

§t u{x > t] + a i u{x ' f) - f (1 - A) & u{x > *) + ••• = o- ( 9 ) 

Truncating Equation ([9| yields the advection-diffusion equation 

^U(x,t) + a-^U(x,t)-v^U(x,t)=0 (10) 
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where v = ah ^ x \ For discontinuous initial data (Hi, U(x, 0) = for x < and 
U(x, 0) = it/{ for x > 0. The analytic solution to (fToTfor i > is then found to be 



erf 



at 



>Avt 



(11) 



where erf(£) is the error function 



erf(C) = - 



2 I"- 
* Jo 



For additional details on this derivation refer to 

The analysis to follow assumes the use of the L\ norm and sets z = x ~ atf, 
S± = ^JAvitf and 82 = ^/Au^tf- Furthermore, assume hi > h,2- As in [11], assume 
that the solution to the modified equation is an accurate approximation to the 
numerical solution so Uh(x, t) = U(x,t). Following a similar line of reasoning as in 
Section [2] gives 



\u hl (x) - u h2 (x)\ 



2 2 \oi 



f 



Si 



Wr - U L \ 

7* 



/11 — \ h 



dz 



Therefore, under the assumption that the three numerical approximations have 

been obtained using the same CFL A, the factor of ^/ at /^ jl \u>r — Ul\ will appear 
in both the numerator and denominator when the ratio of the norms of differences 
is taken. As a result 



||«hi(aO - u h3 {x)\ 




\fh~i — \fh3\ 


\\uh 2 {x) -u h3 {x)\ 




- \/hz\ 



and it is easy to verify that all three approaches to Richardson extrapolation will 
yield convergence at the expected rate of 0.5. 

Remark: Although the results presented in Table [2] use simulations with uniform 
refinement, this is not critical in the analysis for this case. In fact, it is primarily 
the monotone nature of the similarity solution which is responsible for the robust 
nature of the estimates. Uniform refinement was used in order to match the analysis 
for high-order schemes below, where uniform refinement is important. 
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4.2 Arbitrary Order Linear Scheme 

In general, robust results like those represented in Table [5] are not expected. In 
this section we analyze why this is the case and determine a particular strategy 
which yields accurate estimates even in the presence of discontinuities (or other 
self-similar features). We restrict our attention to noncompressive stable p order 
schemes for advection with modified equations of the form 

B 8 B p+1 

Of^t) + OfrVM ~ VhQ^+TU( X ,t) + • • • = (12) 

where r/ h = fjh p and fj is a constant depending on the CFL A. Following the 
analysis in [11] , a simple change of variables is performed to translate into a frame 
of reference traveling with the wave 



z = x — at 

T = t. 



After dropping the higher order terms, ( 12 I becomes 

d -U{ Z ,T) 



QP+ 



dzP+ 



T U(z,r)=0 



(13) 



where is either plus or minus r/^ depending on the value of p. Similarity solutions 
can be sought with similarity variable 



6, 



(14) 



We again assume that the solution of the modified equation is an accurate represen- 
tation of the numerical approximation and set = U . For jump initial condition 
([T]) the solution can then be written in the general form 



S 1 7, b(£h) 



(15) 



where S is an approximation to the jump from —1 to 1 (similar to an error function 
but perhaps with more complex behavior). In general, S will take the form of 
generalized hypergeometric functions which oscillate on one or both sides of the 
discontinuity. In a Richardson style error estimate, norms of the difference between 
two solutions will be used, and so we write 



\u hl (x) - u h2 {x)\ 



o 1 o 6 J 



u L +u R u R - U L 



\Ur - u L \ 
2 



f-OO 


s ( 


/ — 00 





dz. 



Making the change of variables to 
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gives 

/oo 
-co 

Now take the ratio of such norms to arrive at 

W^tn ( x ) - u h 2 ( x ) 
\\ u h 2 (?) - u ha (x) 

In general therefore, the Richardson estimate will depend on the ratio of integrals 
of scaled similarity functions 



rOO 
J — CO 


s(x)-s 


[x p+ < 


l fh7\ 

\j h 2 ) 


dx 


rOO 
J — OO 


s(x)-s 


[x p \ 


sj h 3 ) 


dx 



(19) 



The function S can be an extremely complex object, and so computing this ratio 
in closed form is in general impractical. In fact, the definition of S is often found 
to make this ratio difficult to even estimate numerically due to ill-conditioning 
and finite precision arithmetic. However, for the case of uniform refinement when 
/13 = rh.2 = r 2 hi, this ratio is simply unity. Therefore, for this special case, the 
estimated convergence rate will be given by TL{u hl ,u rhl ,u r 2 hi ) = p/(p+ 1). This 
is the expected convergence rate as discussed in [11] . Notice that such cancelation 
requires both uniform refinement, and that the estimate be performed using dif- 
ferences of successive refinement. Other choices will in general not yield the rate 
p/(p+l). 

To be clear, our analysis shows that for Richardson extrapolation of the form 
H(uh 1 ,u r f ll ,u r 2f ll ) where r is a uniform refinement rate, the a-priori convergence 
rate of p/(p + 1) will be obtained even in the presence of discontinuities. Note here 
that r is simply a uniform refinement, but its exact value is not important in the 
analysis. Other realizations of Richardson extrapolation may not yield this result, 
and the exact value of the computed estimate can vary widely from the expected 
rate. 



s(x)-s\x 



d X - (17) 



rOO 
•J — OO 


s( x )-s 


[x p \ 




dX 


poo 
J — CO 


s(x)-s 


[x p \ 


\J h 3 ) 


dX 



(18) 



4.3 Second Order Linear Scheme 

In order to demonstrate the implications of the analysis in Section |4.2| consider 
the linear second-order upwind method 



n+l n -v 



(^ u i + j(! ~ X )( u i+1 ~ «?-l)J ~ + i(l - A)(u™ - «"_ 2 )j 

(20) 

This is simply a second order unlimited Godunov method. As before, one can 
perform an estimate of the convergence rate using Richardson extrapolation. Again 
we set a = 1, choose a computational domain with [a;^,!^] = [— tt,tv], integrate to 
a time tf =2, and use A = 0.6. The series of approximations is generated using a 
uniform refinement process with a ratio < r < 1 (i.e. hi = rh\ and /13 = r 2 h±), 
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r 


7£ (n hl ,u r/ll ,n r 2 hl ) 


U [u hl ,u r 2 hl ,u rhl ^j 


U (u rhl ,u hl ,u r 2 hl ^j 


l 

2 


0.67 


0.14 


1.63 


2 
5 


0.67 


0.22 


1.73 


1 

3 


0.67 


0.34 


1.60 


2 
7 


0.67 


0.41 


1.49 


1 

4 


0.67 


0.47 


1.35 



Table 3 Estimated convergence rates for the second-order upwind scheme for a solution with 
a discontinuity. The base resolution uses 51201 points, and uniform refinement is carried out 
using a ratio of r. Results using the three independent variants of Richardson extrapolation 
are presented in the various columns. 



and starts with 51201 points in the domain (i.e. hi = 51 2 ^) ). Table [3] shows the 
results using the three basic variants of Richardson extrapolation, and various 
choices of the uniform refinement ratio r. The table shows that the expected rate 
of 2/3 is obtained as advertised when using the appropriate approach. Also shown 
in the table are the type of results that can be experienced if other approaches are 
used. 



4- 3.1 A closer look at this case 

In order to more clearly explain what is going on, we present additional details for 
this case. As discussed above in Section [4. 2| the crux of the matter centers around 
the similarity solution S. For this second-order scheme, the solution can be found 

as 



1 e„ (r (§)) 2 iF 2 (§; |, I; |) - 4. lF2 (I; §, §; |)) 



where r is the Euler Gamma function, and 1F2 is a generalized hypergeometric 
function. Figure[2]shows similarity solutions in the reference frame moving with the 
discontinuity at three resolutions. The grid spacing is essentially a parameter, and 
so we have chosen a normalization hi = 1. Solutions with two grid doublings are 
also shown. Following the analysis in Section |4.2| differences of the three solutions 
will be taken. Figure [3] shows the three sets of differences which are produced for 
the three variants of the Richardson extrapolation error estimate. All three plots 
show the very complex character of the function whose absolute integral is taken. 
The key observation of this paper is presented in Figure[4]where the spatial variable 



is scaled to the common reference variable Xi as suggested in Equation (16 1. For 



the case of uniform refinement when the differences are made as suggested, the 



integrals in the numerator and denominator of ( 19 1 are identical, and the estimate 
follows. In this case the method essentially avoids the need to calculate the actual 
integral and relies on the fact that the ratio is known a-priori as one for any 
similarity function. 
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Fig. 2 Similarity jumps for the unlimited second order scheme in the frame of reference moving 
with the discontinuity. 



1 

0.5 

-0.5 

-1 











S N" S 2N 

S 2N" S 4N 



1 

0.5 



-0.5 











S N" S 4N 

S 2N~ S 4N 



-15 



-10 



-15 



-10 




Fig. 3 Differences of similarity jumps in the frame of reference moving with the discontinuity. 



4.4 Additional demonstration of the theory in ID 



In order to further demonstrate the validity of the theory we have just described, 
we perform a series of tests for linear schemes of increasing order, as well as a 
high-resolution nonlinear TVD scheme. The linear schemes we investigate here 
are upwind biased single step schemes (i.e. in advancing from t n to t n+1 they use 
data from t n only) that are high-order accurate in space and time. Derivation 
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of these schemes follows standard procedures using a Cauchy-Kowalewski pro- 
cess, sometimes called the Lax-Wendroff procedure |16] . Further details on these 
derivations can be found in |17lll8lTl3] . The nonlinear TVD discretization is of the 
high-resolution Godunov type described in [191120] . 

4-4-1 Fourth-Order Linear Scheme 

Consider the linear fourth-order upwind method 

£ Ci%uf +S (22) 

s=-3 

where is a vector of stencil coefficients given by 

5 - 8A 2 + 3A 3 
-37 - 6A + 52A 2 - 9A 3 
146 + 96A - 104A 2 + 6A 3 
-50 - 180A + 80A 2 + 6A 3 ' 
-71 + 96A-16A 2 -9A 3 
7 - 6A - 4A 2 + 3A 3 



(4) = 

144 



Table [4] shows results using the three basic variants of Richardson extrapolation, 
and various choices of the uniform refinement ratio r using ( 22 1 . The test problem 
is again defined a = 1, [a^x^] = [— 7r,7r], tj = 2, and A = 0.6. The first column 
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shows very good agreement between the estimated convergence rate and the ex- 
pected rate of = 0.8. The other two columns indicate that results which are 
difficult to interpret can be obtained for other procedures. Note that the estimated 
rates are slightly higher than the theoretical prediction. This behavior has been 
observed in other studies, and is the result of very slow convergence of high-order 
numerical approximations to the limiting similarity solutions. Despite the large 
number of grid points used in this study, we attribute the slight overestimation of 
the convergence rates in the first column to a lack of sufficient grid resolution. 



r 




Tl ( u h 1 ,U r 2 h 1 ' u rh 1 ' S J 


U (u rhl ,u hl ,u r 2 hl ^j 


l 

2 


0.86 


0.23 


2.32 


2 
5 


0.83 


0.41 


2.10 


1 

3 


0.83 


0.53 


1.91 


2 
7 


0.84 


0.73 


1.60 


1 

4 


0.85 


0.65 


1.34 



Table 4 Estimated convergence rates for the fourth-order upwind scheme for a solution with 
a discontinuity. The base resolution uses 51201 points, and uniform refinement is carried out 
using a ratio of r. Results using the three independent variants of Richardson extrapolation 
are presented in the various columns. 



4-4-2 Sixth-Order Linear Scheme 

Consider the linear sixth-order upwind method 



n+l 



(23) 



where C*- 6 - 1 is a vector of stencil coefficients given by 



C 



(6) 



4320 



-31 + 43A^ 
289 + 24A - 391A 2 - 



- 15A 4 + 3A 5 
30A 3 + 123A 4 



15A & 



-1299 - 324A + 1623A 2 + 360A 3 - 387A 4 
4325 + 3240A - 2675A 2 - 1170A 3 + 615A 4 



-1085 - 5880A + 1505A 2 + 1680A 3 
-2589 + 3240A + 267A 2 
431 - 324A 



+ 27A & 
- 15A 5 
525A 4 - 15A 5 
1170A 3 + 225A 4 + 27A 5 



419A 2 + 360A 3 



33A 4 - 15A & 



30A 3 - 3A 



4 + 3A 5 



-41 + 24A + 47A^ 

Table [5] shows results using the three basic variants of Richardson extrapolation, 
and various choices of the uniform refinement ratio r using ( 23 1 . The test problem 
is again defined using a = 1, = [— 7r,7r], tj = 2, and A = 0.6. The first 

column again shows remarkable agreement between the estimated convergence 
rate and the expected rate of = 0.86. The other two columns indicate that 
results which are difficult to interpret can be obtained for other procedures. As 
in Section |4.4.1| we attribute the slight overestimation of the convergence rates in 
the first column to a lack of sufficient grid resolution. 
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r 


7£ (n hl ,u r/ll ,n r 2 hl ) 


U [u hl ,u r 2 hl ,u rhl ^j 


U (u rhl ,u hl ,u r 2 hl ^j 


l 

2 


0.90 


0.16 


2.95 


2 
5 


0.90 


0.47 


2.35 


1 

3 


0.88 


0.60 


1.93 


2 
7 


0.88 


0.78 


1.20 


1 

4 


0.90 


0.83 


1.24 



Table 5 Estimated convergence rates for the sixth-order upwind scheme for a solution with 
a discontinuity. The base resolution uses 51201 points, and uniform refinement is carried out 
using a ratio of r. Results using the three independent variants of Richardson extrapolation 
are presented in the various columns. 



4-4-3 Second Order Nonlinear Scheme 



Finally, we consider a high-resolution TVD limited scheme. The scheme is a second- 
order MUSCL type scheme using a MinMod limiter applied to the slopes. The 
scheme can be written 



A 



-1 + ^(1 -A)/J 



(24) 



where 



and 



a = MinMod(it" +1 — u™,w™ — u i—x)% 
13 = MinMod(?i" - u - l _ 1 , u ™_ i - w™_ 2 ), 



MinMod(6, c) 



b if \b\ < |c| and bo 
c if |6| > |c| and bo 
if be < 0. 



Note that this is nothing more than a high-resolution Godunov method (see [191 
[20] for details). Table [6] shows the results for the three variants of Richardson ex- 
trapolation convergence estimation. The results in [TT] established that although 



this scheme is nonlinear, the modified equation has solutions of the form (15 1. 
Therefore, our analysis still applies to this case, and the method is expected 
to yield accurate estimates of the convergence rates for uniform refinement case 
1Z (11^ , Urhi i u r 2 h 1 ) • The results in Table[6]show this to be largely true with notable 
qualifications. Notice that the estimates become more accurate with decreasing r. 
As r becomes smaller, the finest resolution becomes relatively more accurate, and 
one might expect an extrapolation estimate to yield results which are more similar 
to those found when comparing computed results to the exact solution. Indeed, as 
seen from the prior results in Tables |3j [4] and [5j even poorly constructed extrapo- 
lation techniques tend to yield somewhat more accurate results as r decreases. In 
addition, there is some significance to the fact that the numerical values in table [6] 
are somewhat less accurate than the results from the second-order linear scheme 
in table [3] The root cause of this observation is not entirely understood, but our 
conjecture is that the discontinuities in the higher derivatives of the similarity 
solution for the MinMod scheme [11] result in slow convergence. Finally, we note 
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r 


7£ (u hl ,u r/ll ,« r 2 hl ) 


U [u hl ,u r 2 hl ,u rhl ^j 


U (u rhl ,u hl ,u r 2 hl ^j 


l 

2 


0.48 


0.48 


0.48 


2 
5 


0.55 


0.55 


0.56 


1 

3 


0.57 


0.57 


0.57 


2 
7 


0.57 


0.57 


0.57 


1 

4 


0.60 


0.59 


0.60 



Table 6 Estimated convergence rates for the high-resolution TVD scheme for a solution with 
a discontinuity. The base resolution uses 51201 points, and uniform refinement is carried out 
using a ratio of r. Results using the three independent variants of Richardson extrapolation 
are presented in the various columns. 



that all the estimates in table |6j including those given by lZ(u hl ,u r 2 hi ,u rhl ) and 
'R-{ u rh 1 , u h 1 > u r 2 h 1 )i appear equally accurate. This is not in contradiction with the 
theory which says that the estimate lZ(ui ll ,u r f ll ,u r 2f ll ) will be accurate but not 
that the others will be inaccurate. In fact this is a similar result to those results for 
the first-order scheme in Table [2] As was the case for the first-order discretization, 
this fortuitous behavior can be traced to the monotonicity of the approximations. 

5 Revisiting the top-hat problem 

Having demonstrated the behavior of Richardson extrapolation error estimation 
for linear jumps in ID, we now return to the introductory example of Section [3j 
Because the second-order wave equation can be written as a system of first-order 
hyperbolic equations, the analysis presented in the previous sections is expected 
to have applicability here also. We therefore conduct a similar set of studies to 
those presented in Section [4] We perform extrapolation estimation to a series of 
computations that use a uniform refinement process and present results for the 
three variants of Richardson extrapolation while varying the refinement ratio r. 
Tables [7| through [10| present these results. 



r 


K (u hl ,u rhl ,u r 2 hl ) 




K (u rhl ,u hl ,u r 2 hl ^ 


l 

2 


0.44 


0.40 


0.50 


2 
5 


0.47 


0.43 


0.53 


1 

3 


0.47 


0.44 


0.53 


2 
7 


0.47 


0.44 


0.54 


1 

4 


0.47 


0.44 


0.53 



Table 7 Estimated convergence rates for the first-order upwind scheme and the top-hat prob- 
lem. The base resolution has 400 grid points per dimension and uniform refinement with ratio 
r is carried out. Results using the three independent variants of Richardson extrapolation are 
presented in the various columns. 
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r 


TZ [u hl ,u rhl ,u r 2 hl ^j 


TZ (u hl ,u r 2 hl ,u rh ^ 


TZ (u rhl ,u hl ,u r 2 hl ^j 


l 

2 


0.53 


-0.27 


1.92 


2 
5 


0.56 


-0.08 


2.02 


1 

3 


0.55 


0.06 


1.90 


2 
7 


0.57 


0.17 


1.87 


1 

4 


0.57 


0.25 


1.73 



Table 8 Estimated convergence rates for the second-order upwind scheme and the top-hat 
problem. The base resolution has 400 grid points per dimension and uniform refinement with 
ratio r is carried out. Results using the three independent variants of Richardson extrapolation 
are presented in the various columns. 



r 


TZ (u hl ,u rhl ,u r 2 hl ^ 


TZ (u hl ,u r 2 hl ,u rhl ^ 


TZ (u rhl ,u hl ,u r 2 hl ^ 


l 

2 


0.62 


0.52 


0.79 


2 
5 


0.63 


0.55 


0.78 


1 

3 


0.63 


0.57 


0.75 


2 
7 


0.63 


0.58 


0.76 


1 

4 


0.63 


0.58 


0.74 



Table 9 Estimated convergence rates for the nonlinear high-resolution upwind scheme and 
the top-hat problem. The base resolution has 400 grid points per dimension and uniform re- 
finement with ratio r is carried out. Results using the three independent variants of Richardson 
extrapolation are presented in the various columns. 



r 


TZ (u hl ,u rhl ,u r 2 h \ 


TZ (u hl ,u r 2 hl ,u rhl ^ 


TZ (u rhl ,u hl ,u r 2 hl 'j 


l 

2 


0.64 


-0.50 


3.41 


2 
5 


0.65 


-0.11 


3.36 


1 

3 


0.65 


0.11 


2.85 


2 
7 


0.66 


0.25 


2.63 


1 

4 


0.66 


0.34 


2.33 



Table 10 Estimated convergence rates for the fourth-order upwind scheme and the top-hat 
problem. The base resolution has 400 grid points per dimension and uniform refinement with 
ratio r is carried out. Results using the three independent variants of Richardson extrapolation 
are presented in the various columns. 



The results show that the conclusions derived from the simple ID analysis hold 
even for this more complex situation. In particular, the estimated rate of conver- 
gence found when using uniform refinement and using 7?. (u^ , u,,/^ , u^/^ ) yields 
results that are in reasonable agreement with the a-priori expected rate p/(p+ 1). 
However, the two other variants are found to yield wildly varying results for both 
the second-order and fourth-order cases. As for the ID case, the first-order and 
high-resolution methods yield monotone (or nearly monotone) approximations, 
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which leads to fortuitous cancelation. As a result the estimated convergence rates 
are nearly correct for all three methods. This behavior is not a general result, and 
should not be relied upon in practice. One additional interesting point is that as the 
separation between the grid spacing of the three numerical solutions increases, the 
estimated convergence rate seems to converge toward the a-priori expectation. For 
example, the estimates in the second and third columns in Table [To| are woefully 
inaccurate, but in both cases the accuracy is improving as the relative separation 
increases. This may be one reason why practical studies, such as those in 8 , have 
used Richardson extrapolation with extremely refined final simulations. 



6 Conclusions 

We have provided an in depth investigation of Richardson extrapolation error es- 
timation for linear hyperbolic wave propagation in the presence of discontinuous 
solutions. The discussion is motivated by an example of 2D scalar wave propaga- 
tion that has been taken from the literature. Numerical methods of orders 1,2,4, 
and a high-resolution (nominally 2 n< ^ order) variant are used to produce approxi- 
mations at three resolutions. Richardson extrapolation error estimates found using 
these approximations are found to produce results that can vary widely from the 
expected convergence rate of p/(p + 1) where p is the order of the method. One ex- 
ception to this wide variation is observed if the estimate is determined in a specific 
manner. 

In order to understand this behavior, a simpler model problem of ID advection 
is posed as an appropriate surrogate, and a detailed analysis is presented. The 
analysis uses the solution to the modified equation to elucidate the difficulty found 
in practice. In addition, the analysis reveals that one particular realization of 
the technique reproduces the a-priori convergence rates even in the presence of 
a discontinuity or other similarity type behavior. The key elements are shown 
to be the use of uniform refinement, and that the comparisons inherent in the 
Richardson estimator are performed in one specific manner as 1Z(uf ll ,u r f ll ,u r 2f ll ) 
where r is the uniform refinement ratio. This result was then demonstrated in ID 
for a number of discretizations ranging in order from first-order to sixth-order, and 
for the motivating 2D discretizations ranging from first-order to fourth-order. In 
addition, results were presented for nonlinear high-resolution schemes. The results 
were found to be in good agreement with the theory. 

The analysis presented in this work is an interesting result that shows why 
estimated convergence rates from Richardson extrapolation for problems with dis- 
continuities have been found to be difficult to interpret. In addition, it shows that 
Richardson extrapolation can be used to obtain predictable a-priori convergence 
rates even for simulations of systems with linear discontinuities provided two key 
elements are satisfied. First, uniform refinement must be used. Second, the esti- 
mate must be obtained as TZ{uj ll ,u r } ll ,u r 2f ll ) in the notation of this paper. The 
result is in fact more general than discontinuities and includes other self-similar 
features such as corners. Future work will include investigating the techniques dis- 
cussed here for nonlinear equations such as the Euler equations, as well parabolic 
problems or dispersive wave propagation problems. 
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